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We study the dynamic properties of a model for wetting with two competing adsorbates on a planar 
substrate. The two species of particles have identical properties and repel each other. Starting with 
a flat interface one observes the formation of homogeneous droplets of the respective type separated 
by nonwet regions where the interface remains pinned. The wet phase is characterized by slow 
coarsening of competing droplets. Moreover, in 2+1 dimensions an additional line of continuous 
phase transition emerges in the bound phase, which separates an unordered phase from an ordered 
one. The symmetry under interchange of the particle types is spontaneously broken in this region 
and finite systems exhibit two metastable states, each dominated by one of the species. The critical 
properties of this transition are analyzed by numeric simulations. 



INTRODUCTION 



Wetting phenomena are observed in a large variety of 
situations where an inert surface is exposed to a bulk 
phase such as a gas or a liquid. Depending on external 
parameters like chemical potential and temperature the 
internal forces between the particles and the solid may 
lead to the formation of a thin layer of a different thermo- 
dynamical phase, a so-called wetting layer The mor- 
phology and the thickness of the layer depends on how 
the free-energy contributions at the interfaces between 
the solid, wetting layer and the bulk phase balance one 
another. If the bulk phase is thermodynamically favor- 
able the wetting layer remains bound to the surface and 
is characterized by microscopically finite average width. 
However, approaching the point where the gas phase and 
the wetting layer coexist in the bulk, the system under- 
goes a wetting transition. Beyond this transition the wet 
phase becomes more favorable in the bulk so that the 
layer grows, eventually reaching a macroscopic size. 

The phase diagram of a system with wetting layers 
could be rather complex exhibiting a variety of surface 
phase transitions,prewetting phenomena, and multicrit- 
ical behavior 0, Hf. For example, if the temperature is 
varied while moving along the coexistence curve of wet- 
ting layer and bulk phase, a transition may take place at a 
temperature Tw beyond which the thickness of the layer 
becomes infinite. This transition, known as continuous 
wetting, is usually first order, although in certain models 
the transition is continuous. On the other hand, varying 
the chemical potential difference between the two phases 
and moving towards the coexistence curve at T > T\y, a 
different type of transition takes place, which is referred 
to as complete wetting. 

In many situations it is reasonable to assume that a 
statistically stationary layer is in thermal equilibrium 
with its environment. In the 1980's many authors studied 
such wetting transitions experimentally and theoretically 
within the framework of equilibrium statistical mechan- 
ics (for a review, see 0). Within this approach, a wet- 
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FIG. 1: (color online) Sketch of a wetting process with com- 
peting adsorbates: Two types of particles A and B, which 
strongly repell each other, are adsorbed on top of an inert sub- 
strate, forming a wetting layer of mutually avoiding droplets. 



ting transition is usually modeled as the unbinding of an 
interface from a wall. The interface configuration is de- 
scribed by a height function /i(x) above the point x on 
the substrate. The model is then defined in terms of an 
effective energy functional [1] 



d d x 



: (Vh) 2 + V(h(x)) 



(1) 



where a is the effective surface tension of the interface, 

V (h) is a potential accounting for the interaction between 
the wall and the interface, and d is the interface dimen- 
sion (usually d = 2). In the nonwet phase the potential 

V contains an attractive component which binds the in- 
terface to the wall. Assuming thermal equilibrium, the 
probability of finding the interface in a certain configu- 
ration is then given by the canonical distribution 



P[h]~exp(-f3£[h}). 



(2) 



As the parameters describing the system are varied, 
the attractive component of the potential may become 
weaker so that it is no longer able to bind the interface, 
leading to a wetting transition. 

If one is interested not only in static properties but also 
in time-dependent features such as dynamical roughen- 
ing, one usually introduces a stochastic Langevin equa- 
tion that reproduces the equilibrium distribution ([2|) in 
the limit t — ► oo. Assuming short-range interactions and 
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keeping only the most relevant terms in the renormal- 
ization group sense, one is led to the Edwards- Wilkinson 
equation with a potential [j| 



dh(x,t) 
dt 



, , . dV(h(x,t)) , s 



where £(x, i) is a zero-average Gaussian noise field with 
a variance 



(C(x,t)c(x',t')) = 2r^- 1 (x- 

and a noise amplitude T = ksT. 



x')6(t-t'), (4) 



In contrast to the canonical ensemble, the definition 
of the model as a dynamical process allows one to go 
beyond the realm of equilibrium thermodynamics and to 
study the influence of nonequilibrium effects, which man- 
ifest themselves as a violation of detailed balance. In 
fact, during the past decade the study of wetting transi- 
tions far from equilibrium e merg ed a s a new sub-topic 
of the field [1, 0,1, S, EE El El El El- In these theo- 
retical studies it was shown that nonequilibrium wetting 
transitions may differ significantly from their equilibrium 
counterparts, exhibiting e.g. different types of critical be- 
havior and new macroscopically observable physical phe- 
nomena. 

In the present work we investigate what happens if the 
surface is exposed to a gas consisting of two different 
types of particles, say A and B. Similar situations, but 
with a focus on catalytic processes, were investigated re- 
cently in The two species of particles are assumed 
to repel each other strongly, leading to the formation of 
competing droplets of either type, as sketched in Fig. Q] 
In our model the repelling force is implemented as a dy- 
namical constraint that deposition is prohibited if it were 
to cause a direct contact between an A and a B particle. 
Moreover, the two types of particles are assumed to have 
identical physical properties, establishing ai<-*B sym- 
metry in the model. Certainly it is difficult if not impos- 
sible to realize these assumptions experimentally. How- 
ever, here we are primarily interested in the theoretical 
question of how such a symmetry influences the universal 
critical behavior of wetting transitions. For simplicity the 
present work is restricted to the case of detailed balance, 
where methods of equilibrium statistical mechanics can 
be applied. However, the generalization to the nonequi- 
librium case is straightforward. 

The main results of this paper are the following: 

(i) The competition of A and B islands slows down 
the dynamics significantly, particularly in the unbound 
phase, where an extremely slow coarsening process is ob- 
served. 

(ii) In the case of a one-dimensional substrate the phase 
structure is the same as in the single-species case. 

(iii) In two (and higher) dimensions an additional tran- 
sition line divides the bound phase into two parts. 
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FIG. 2: (color online) Dynamics of the two-species model. 
Particles are deposited on the substrate at rate qo (a) and on 
islands of the same type at rate q [(b) and (g)]. Moreover, 
they evaporate from the edges of terraces at rate r (c) and 
from the middle of plateaus at rate p (h). Particles of different 
types are not allowed to touch each other [(d) - (f)]. 



(iv) For a growth rate q = this additional transition 
is found to belong to the same universality class as the 
kinetic Ising model. For q = 1 one observes an unusual 
type of transition. 



DEFINITION OF THE MODEL 



Single-species model 



Before introducing the two-species model let us briefly 
recall the single-species model introduced in [||. The 
single-species model is defined as a solid-on-solid growth 
process on a flat substrate represented by a ci-dimensional 
hypercubic lattice with N — L d sites. The configuration 
of the wetting layer is described by an interface without 
overhangs, meaning that each lattice site i of the sub- 
strate is associated with an integer height variable hi. 
The interface evolves by deposition and evaporation of 
particles restricted by two constraints. On the one hand 
the interface obeys the so-called restricted solid- on- solid 
(RSOS) condition that the heights of neighboring lattice 
sites i,j may differ by at most one unit: 



hA < 1. 



(5) 



This constraint was first introduced in [16[ and imposes 
an effective surface tension. On the other hand, the 
chemically inert substrate is modeled as a hard-core wall 
by imposing the dynamical constraint 



hi > 0. 



(6) 



The model evolves by random-sequential dynamics, i.e., 
a site of the lattice is randomly selected and one of the 
following processes is carried out, provided that it does 
not violate the constraints and ©: (i) deposition on 
the substrate at rate qo , (ii) deposition on top of islands 
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at rate q, (iii) evaporation from the middle of plateaus at 
ratep, (iv) evaporation at the edges of plateaus at rate r. 

If the resulting configuration were to violate the con- 
straints §5§ or ^ the attempted move is rejected. As 
usual, the time scale is fixed by choosing one of the rates, 
e.g., r — 1. For p ^ 1 the process defined above can be 
shown to violate detailed balance in the stationary state, 
while for p — 1 and q < I the stationary state is given by 
a Boltzmann-Gibbs distribution. 

The wetting transition in this model can be understood 
as follows. Far away from the wall a free interface prop- 
agates with a certain average velocity 4t{h) = v(q,p) 
which depends on the growth rate q. Varying the growth 
rate this velocity changes sign at a well-defined thresh- 
old q c {p)- For v(q,p) < the interface moves backwards 
until it reaches the substrate where it continues to fluc- 
tuate in a stationary bound state, while for v(q,p) > 
the interface detaches from the bottom layer. 



Two-species model 

In the present work we generalize the model introduced 
in Q as to describe deposition and evaporation of two 
different types of particles, labeled by A and B. Both 
species of particles are completely symmetric and obey 
the same dynamic rules as in the original model. In ad- 
dition we require that the two species repel each other, 
implemented by the dynamical constraint that particles 
of different types are forbidden to be in contact with each 
other. For example, a particle of type A must not be de- 
posited on top of a B and vice versa. Similarly, particles 
of different types cannot be deposited at adjacent sites 
(see Fig. [2]). Starting with a flat interface at zero height, 
these dynamical rules ensure that at each lattice site the 
wetting layer consists of only one type of particles, lead- 
ing to homogeneous droplets, consisting either of A or B 
particles. Note that this model is completely symmetric 
under an interchange A «-> B, establishing a global Zi 
symmetry. 

The present study is restricted to the case p = 1, where 
stationary bound states obey detailed balance so that 
methods from equilibrium statistical mechanics can be 
applied. A full analysis of the genuine nonequilibrium 
case p ^ 1 will be presented elsewhere. 

As sketched in Fig. [3j the generalized model can be 
implemented numerically in a very elegant way by repre- 
senting a wetting layer of type A by positive and a layer 
of type B by negative heights. As before, hi = stands 
for an unoccupied site where the substrate is exposed to 
the gas phase. Compliance with the additional rule that 
particles of different types may not be in contact with 
each other is then consistently ensured by extending the 



TABLE I: List of the probabilities for deposition and evapo- 
ration in a local update, normalized by n = max(r + q, 2 go). 
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FIG. 3: (color online) Numerical implementation: Droplets 
consisting of A and B particles are represented by positive 
and negative heights, respectively. 



RSOS constraint to negative heights. Using this rep- 
resentation, the generalized model can be simulated by 
carrying out the following local update rule: 

1. A site i is chosen at random. 

2. A deposition or evaporation process is selected ac- 
cording to the probabilities listed in Table |TJ 

3. The attempted update is rejected if the resulting 
configuration violates the restricted solid-on-solid con- 
dition ([5]), otherwise it is carried out. 

4. The time variable is incremented by At = —j^, 
where n = max(r + <?, 2 go)- 

The successive execution of N local updates is referred 
to as one Monte Carlo step (MCS). As before, the overall 
time scale can be fixed by choosing one of the rates, e.g. 
by setting r = 1. 

The representation of A and B layers by positive and 
negative heights is not only technically useful. As we 
will see below, it is also instructive in analytical consid- 
erations, as it allows one to relate the two-species model 
to the original single-species model. 

In what follows we are primarily interested in the phys- 
ically relevant case of a two-dimensional substrate be- 
cause this is the lowest dimension where the competition 
of the two particle species leads to an additional phase 
transition in the pinned region. In all numerical simu- 
lations we use periodic boundary conditions in order to 
minimize finite size effects. 



PHASE DIAGRAM 

To describe and analyze the model, we mainly consider 
three different order parameters. The first quantity is the 
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FIG. 4: (color online) Phase diagram of the two-species model 
for p = 1 in 2+1 dimensions, comprising a symmetric phase 
(A), a symmetry-broken phase (B), and a rough phase (C/C). 

occupation balance, defined by 

1 N 

t=i 

Since particles of type A are represented by positive and 
particles of type B by negative occupation numbers hi, 
the occupation balance is a measure of the surplus of 
one type of particles and thus indicates a broken A <-> B 
symmetry. 

The second quantity is the density of unoccupied sites 
aK°) 1 N 

i=l 

This parameter indicates how strongly the interface is 
pinned to the bottom layer. In particular, for large de- 
position rates, where the surface is covered by compet- 
ing islands which are separated by lines of unoccupied 
sites, this order parameter is a measure of the length of 
the perimeter of the islands and thus allows us to draw 
conclusions about their size and the roughness of their 
borders. 

Finally, we are also interested in the interface width 

« = J^En 2 -(^En) 2 ( 9 ) 

y i i 

which quantifies the roughness of the interface. These or- 
der parameters allow one to identify the following phases 
(see Fig. HJ: 

(A) Symmetric phase: For q < 1 and small values of go 
the interface fluctuates close to the wall, forming small 



short-lived islands of either type. The average occupation 
balance is zero and the width of the interface saturates 
at a finite value. 

(B) Symmetry-broken phase: For large values of go 
and q < 1 one observes coarsening patterns of com- 
peting islands. The islands have a pancake like shape 
and the interface width stays finite. Finite systems have 
two metastable states, that each feature a surplus of one 
of the particle types and become stable in the thermo- 
dynamic limit. The model has a stationary symmetry- 
broken state with a nonzero occupation balance. This 
phase does not exist in 1+1 dimensions. 

(C,C) Rough phase: As free parts of the interface tend 
to grow for q > 1, this phase is characterized by the com- 
petition of islands of the two particle types. As these 
islands coarsen, one observes a monotonic increase of the 
average interface width. This coarsening process turns 
out to be extremely slow. 

In finite systems a state with only one type of particles 
is reached after a characteristic time scale. In the limit of 
large systems this time scale almost certainly grows ex- 
ponentially or super-exponentially with the system size. 
After the expulsion of one of the particle types the inter- 
face detaches from the bottom layer, entering a moving 
non-stationary state. Note that this behavior differs sig- 
nificantly from the one of the single-species model, where 
the interface detaches immediately. 

The phases (A) and (B) are separated by a phase tran- 
sition at a well-defined critical threshold. This raises the 
question whether this line extends to the upper part of 
the phase diagram, separating the rough phase into two 
parts (C) and (C). As will be discussed below, these re- 
gions are quite different in character. However, we were 
not able to locate a transition line between (C) and (C) 
in a robust and reproducible way. Therefore we believe 
that the apparent phase boundary for q > 1 is a precur- 
sor of the (A)-(B) transition rather than a line of genuine 
phase transition. 



CLASSIFICATION OF THE PHASES BY MEANS 
OF THE BINDER CUMULANT 

In phase (B) one observes coarsening and pancakelike 
islands which compete one another. Finite systems have 
two metastable states, which are each dominated by one 
of the particle types. This leads to a spontaneous break- 
ing of the A <-> B symmetry, in a similar way as in fer- 
romagnets. In finite systems the symmetry-broken state 
is only metastable against fluctuations. Therefore one 
observes occasional flips from an A-dominated to a B- 
dominated state and back. Fig. [5] shows a sequence of 
six snapshots of such a flip extending over 4000 Monte 
Carlo updates. As can be seen, a flip occurs whenever 



FIG. 5: (color online) Six snapshots of a 50 x 50 Lattice during a flipping process with qo = 1.17 and q = 0.82 from a state 
dominated by red (bright) to a state dominated by blue (dark) particles. The brightness indicates the height, white stands for 
unoccupied sites. 
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FIG. 6: (color online) Occupation balance b and density of un- 
occupied sites in a system with 50 x 50 sites for q — 1.17 
and q = 0.82. Every data point was averaged over 100 MCSs. 
The system flips between two metastable states. Obviously 
the flipping is correlated with a small increase of the density 
of unoccupied sites. 
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FIG. 7: (color online) The characteristic time interval Tj be- 
tween two flipping processes depends exponentially on the 
system length L. The two metastable states of the symmetry- 
breaking phase become stable in the thermodynamic limit. In 
2+1 dimensions the two-species model exhibits spontaneous 
symmetry-breaking. 
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minority islands generated by fluctuations reach a critical 
size that suffices to displace the majority type. 

The repeated flipping in finite systems can be observed 
by monitoring the occupation balance b. A typical plot 
of the behavior in this phase is shown in Fig. [51 As can 
be seen, the occupation balance b fluctuates around one 
of two values ±6* which characterize the two metastable 
states. This behavior is intermitted by flipping processes 
from one of the states to the other. As shown in Fig. [7] 
the average time between two flipping processes Tf grows 
exponentially with the system size. Thus in the limit 
L — > oo the model has two thermodynamically stable 
ground states that spontaneously break the symmetry 
under interchange of the two particle types. 

Fig. [5] reveals another interesting fact: Obviously, the 
flipping is weakly correlated with a transient increase 
of the density of unoccupied sites. To understand this 
increase qualitatively let us again consider the snap- 
shots shown in Fig. [5] While the system is in one of 
the metastable states the number of unoccupied sites is 
correlated with the number and the size of the minor- 



ity islands, as most unoccupied sites are arranged along 
their perimeter. However, during the flipping the sys- 
tem is dominated by two big islands, one of each particle 
type, that are separated by a fissured, but preferentially 
straight, domain wall, along which the unoccupied sites 
are located. For this, the number of exposed sites (gray 
pixels) is maximal during the flipping process (see Fig. [5j. 

The two phases are associated with different character- 
istic probability distributions of the occupation balance. 
In the symmetric phase (A) a finite system exhibits a 
Gaussian distribution centered around zero, while in the 
symmetry-broken phase (B) the distribution shows two 
peaks localized at ±6*. This qualitative difference can be 
studied by measuring the Binder cumulant 17 1 



3 (b 2 ) 2 



(10) 



In the symmetric phase, where b is normally distributed, 
the cumulant vanishes. In the symmetry-broken phase, 
where the occupation balance fluctuates only slightly 
around ±6*, the Binder cumulant should tend to the 
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FIG. 8: (color online) Binder cumulant U as a function of 
the parameters q and qo in a 2+1-dimensional system with 
TV = 50 2 sites averaged over T = 2x 10 5 MCSs after an equally 
long relaxation time. The same data are shown as a density 
plot in the background of Fig. U Orange (A) corresponds 
to the symmetric phase, violet (B) to the symmetry-broken 
phase and blue (C and C) to the growing phase. 
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(11) 



Interestingly the Binder Cumulant can also be used to 
identify regions in which the interface is propagating at 
a constant velocity v. Such regions can be expected to 
exist in the region q > 1. In this case bt ~ ivt, such 
that one finds 



U 



:ir i/T fi(vt)*dt* 



(12) 



As demonstrated in Fig. [5] these values are nicely re- 
produced in numerical simulations and can be used to 
distinguish between the phases (A) and (B). For q > 1, 
the interface of a finite system detaches from the bottom 
layer after a characteristic time Ta- In region (C) this 
time is much smaller than the simulation time such that 
the Binder cumulant takes the value from Eq. (fF2| . In 
region (C), however, the cumulant remains close to zero. 
Note that this sudden change in the Binder cumulant 
between (C) and (C) does not necessarily mean that a 
phase transition takes place, rather it simply indicates 
that the coarsening process in region (C) is so slow that 
the actual simulation time does not suffice to observe the 
detachment of the interface. 



STATIONARY EQUILIBRIUM STATE 



Detailed balance 



For p = 1 and q < 1 the model is in the bound phase 
and evolves towards a fluctuating stationary state. As 
will be shown below, this stationary state obeys detailed 
balance so that methods of equilibrium statistical me- 
chanics can be applied. More specifically, it turns out 
that the stationary state is characterized by a canonical 
distribution, i.e., the probability of finding an interface 
configuration {hi]f =1 which is consistent with the RSOS 



condition ([5]) is given by 



P({h l }l 1 ) = —e X p[ 

AN 



-w({M<=i)] 

AN 

Here l~L{{hi}f =l ) is an energy functional and 
" -H({hi}? =1 )] 



'N 



exp[ 



(13) 



(14) 



is the partition sum over all configurations {hi}^ obey- 
ing the RSOS constraint. The energy functional is of the 
form 



w({Mf=i) 



N 

£ 

i=l 



V(hi) 



(15) 



meaning that it associates with every height value hi 
a potential energy V(hi) that depends on the rates qo 
and q. Following [7| this potential is given by 



V(h) 



-HqM h = 
-\h\\n(q) h^O 



(16) 



Note that this potential is symmetric under the change of 
sign h — > —h, reflecting the A <-> B symmetry. Inserting 
Eqs. (HSJ) and $16j) into Eq. (JT3J) yields the distribution 

(17) 



P({Mf=x) = ^-ffG^W) (q/q )^ S ^. 
AN 



Following [7J] one can easily show that the distribu- 
tion jig is in fact stationary and obeys detailed balance. 
According to Table Q] the deposition of a particle on the 
substrate or on top of an island changes the probabil- 
ity of a configuration by a factor of qo or q, respectively. 
Likewise evaporation contributes a factor 1/go or l/q. 
As the rates for deposition are go and q and the rate for 
evaporation is 1, this state obeys detailed balance and is 
therefore stationary provided that q < 1. 



Origin of the phase transition 

Regarding Eq. (jTTJ) , the existence of a phase transi- 
tion seems to be surprising as there is no explicit inter- 
action, not even a short ranged one, between the sites 
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of the lattice. In fact, the only mechanism that leads to 
correlations between the sites and therewith to a phase 
transition in this model is the RSOS constraint. 

To understand this mechanism, it is instructive to con- 
sider the simple case q = 0, where the dynamics is re- 
stricted to a single monolayer, i.e. hi S {0,±1}. In this 
case the potential (flBj) is given by 



V(h) 







ln(go) 



h = 



(18) 



and the probability for finding a configuration with hi £ 
{0, ±1} that is consistent with the RSOS constraint (J5]) 
reads 



P({hi} 



N 



l) — 7 % 



1 Ar(+)+Ar<- 

7 % 
Zjv 



— a w 



_AT(0) 



(19) 



Here denotes the total number of sites with hi — 

±1, so that N = AT(+) + AT(-) + The occupation 

balance is then given by b = -±.(AT (+) - #(->). Obvi- 
ously, for small q the system prefers sparsely occupied 
configurations whereas for large q a configurations with 
a high density of islands become more likely. As the 
probability for a specific configuration depends only on 
the number of occupied states + the behavior 
of the system is determined by the number of possible 
configurations m(N^ + \ N^) for given particle numbers 
N^ + > and N^~\ This multiplicity of states plays a key 
role for the phase transition between phases (A) and (B). 

Without the RSOS constraint the multiplicity would, 
regardless of the dimensionality of the lattice, just be 
given by 



m(JV (+) ,iV (_) ) = 



AT(+)!iV(-)!iV( )!' 



(20) 



Therefore, without the RSOS condition, the most prob- 
able configurations would be characterized by s» 
, irrespective of the value of go- 

The influence of the RSOS constraint on the multiplic- 
ities m(N( + \ N^) depends strongly on the dimension- 
ality of the lattice. This is demonstrated in Fig. [51 where 
the multiplicities for systems with N = 25 sites are com- 
pared for the unrestricted and the restricted case in 1+1 
and 2+1 dimensions. As can be seen, the RSOS condi- 
tion reduces the area where m is nonzero since islands of 
different types have to be separated by unoccupied sites 
at zero height. On the other hand, in 2+1 dimensions 
the influence of the RSOS constraint is so strong that 
the maxima are shifted away from the diagonal to the 
edges of the diagram. Therefore, if qo is larger than some 
critical value q$ , the system is preferentially driven into 
one of these configurations at the edges, explaining why 
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FIG. 9: (color online) Multiplicity of states for all possible 
values of JV' + ^ and TV*- - ' of systems with 25 sites. Light Gray 
stands for multiplicity zero. 



the system prefers a symmetry-broken state. Fig. O also 
demonstrates why the symmetry breaking is not observed 
in 1+1 dimensions, where the influence of the RSOS con- 
straint is apparently not strong enough to shift the max- 
ima of the multiplicities to the edges. 

For q < 1 higher layers become accessible but they 
are occupied only rarely so that the islands retain their 
pancakelike shape. Therefore the situation can be ex- 
pected to be qualitatively similar. Nevertheless, particles 
in higher layers tend to stabilize the particles in the layers 
below, increasing the effective growth rate. Therefore the 
critical value q$ should decreases with increasing q. As 
can be seen in the phase diagram (Fig. [4j) , this is indeed 
what the numerical simulations show. 



CRITICAL PROPERTIES 



In this section we study the dynamical critical behavior 
along the transition line between the phases (A) and (B) 
and at its upper terminal point. Here the two order pa- 
rameters defined above require a separate treatment with 
different initial conditions. For the occupation balance bt 
one has to start with a monolayer of only one species of 
particles, say A-particles so that bo — 1. With this initial 
condition the occupation balance b t is expected to decay 
at criticality as 



r 



On the other hand, the density of exposed sites p t 
quires to start with a fiat interface so that p °' = 1. 
criticality this quantity is expected to decay as 



(21) 

re- 
At 



(0) 



Pi 



(o) 



P 



(0) 



(22) 



while the occupation balance vanishes for all t. 



Moving away from the critical line into the symmetry- 
broken phase (B) the occupation balance becomes 
nonzero and flips between the values ±6*. One can there- 
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fore use 



b* = b*(q ,q) = (\b\), 



(23) 



as a magnetizationlike order parameter to characterize 
the transition between the two phases. Suppressing flips 
by choosing sufficiently large system sizes and varying go 
while keeping q < 1 fixed b* is found to increase as 

b*(q ,q)~[qo-q C o(q)f- (24) 

By seeking for power laws and by measuring the order 
parameter b* we first determined various critical points 
gg which are listed in Table HT1 



Transition for q — 



Measuring the dynamic critical exponents a and 5 at 
the lower end of the transition line go — <Zo(0)> Q — wc 
find 



5 = 0.058(4) and 



0.5(1). 



(25) 
,(o) 



The density of exposed sites relaxes vs the value poJ = 
0.384(5). Moving into the ordered phase (B) with q = 
fixed the magnetizationlike order parameter b* obeys 
Eq. (|24|) with a critical exponent of 



= 0.125(5). 



(26) 



This suggests that at this point the transitions belongs 
to the universality class of the two-dimensional kinetic 
Ising model with heat bath dynamics, which is char- 
acterized by the exponents (3 — 1/8, z ~ 2.125, and 
a = (3 1 ' z w 0.0588 [19(. The ordered phase of the kinetic 
Ising model is known to be characterized by coarsening 
domains with an average size growing as 

fV2 18]. One 

would therefore expect that the density of exposed sites 
pf^ , which in the present case is a measure of the density 
of domain walls, decays as t~ 1 / 2 , which is consistent with 
our results. 

This result seems to be reasonable, as at q — 0, where 
the height values are restricted to hi E {0, ±1}, our model 
is very similar to a kinetic Ising model. The behavior is 
governed by the competition of bulk noise, caused by 
evaporation and deposition, and the ordering influence 
of the RSOS constraint. The rate go acts as an effective 
inverse temperature. Similar behavior was observed in 
other binary mixture models on two-dimensional lattices 



Transition line for < q < 1 



For < g < 1 the coarsening islands retain a pan- 
cakelike shape since higher layers are exponentially sup- 
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FIG. 10: (color online) Finite size scaling at q = 0.2 to support 
the hypothesis that the transition is nonuniversal. The dashed 
gray line corresponds to a slope of /3 = 1/8. 



TABLE II: Numerical values for the critical points along the 
phase boundary q^(q) between the unordered phase (A) and 
the symmetry-broken phase (B) and various critical exponents 
a, <5, and /3 along that line obtained from simulations of sys- 
tems with 100 x 100 sites. 



q 





0.2 


0.4 


0.6 


0.8 1 


qo(q) 


2.062(5) 


1.787(5) 


1.560(5) 


1.360(5) 


1.181(5) 1 


a 


0.5(1) 


0.5(2) 


0.5(2) 


0.6(2) 


0.6(2) 


8 


0.058(5) 


0.063(5) 


0.068(5) 


0.074(5) 


0.098(5) - 


13 


0.125(5) 


0.14(1) 


0.16(2) 


0.19(3) 


0.21(3) - 



pressed, as can be seen from the stationary equilibrium 
distribution (|17p . This means, that the third dimension 
(height) is switched off so that the resulting process is 
still effectively two-dimensional. One might therefore ex- 
pect that the transitions along the entire line between the 
phases (A) and (B), except for the upper terminal point, 
belong to the universality class of the two-dimensional 
kinetic Ising model. 

Our simulations suggest that this is not the case, but 
rather that all three exponents vary continuously when 
q is increased. The determination of critical exponents 
is always a delicate issue, especially when, as in the case 
of a and /?, the power laws contain additional unknown 
quantities. Being skeptical of these results we have care- 
fully estimated the error bars. The results of our analysis 
are shown in Table HT1 

The critical points listed in Table [Til were estimated 
in two ways: first by searching for a power law decay of 
the form Eq. (|21j) while varying go for fixed values of g. 
Second by fitting a power law of the form (|24|) to data 
obtained by stepwise increasing go and measuring b* after 
sufficiently long relaxation times. Both methods yield 
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FIG. 11: (color online) Magnetization exponent /3 at q = 0, 
q = 0.2, and q — 0.4 for various system sizes obtained from 
finite size simulations. For system sizes L > 60 the main 
error source is the uncertainty in the position of the phase 
boundary go In- 



consistent results near the lower and upper end of the 
transition line. In the middle, i.e., around q « 0.5, the 
q^{q) values obtained by the second method are slightly 
smaller. We consider the first method to be more reliable 
as it is less affected by finite size effects. 

The exponent S can be determined with high accuracy 
although its value is fairly small. It is the most reliable in- 
dicator for a continuous variation of the exponents along 
the transition line. 

Concerning the exponent a the by far dominating 
source of errors is the determination of the asymptotic 
value of the density of exposed sites, which seems to take 
the value — 0.38(1) along the entire transition line. 
Although we see a systematic dependency of the expo- 
nent a on the rate q this dependency is not significant. 

The most problematic exponent is (3. First, the uncer- 
tainty in the position of the phase boundary q§ (q) limits 
the reachable accuracy. In addition, when q is increased 
the correlation length in phase (B) grows, increasing the 
influence of finite size effects. The feasible system sizes 
are in turn limited by the necessity for very long relax- 
ation times. In order to handle both types of errors we 
performed additional control simulations with stepwise 
decreased deposition rate qg and simulations with differ- 
ent system sizes L < 100. The results of our finite size 
simulations are shown in Fig. \W\ and Fig. QT] The esti- 
mated exponent decreases with increasing system sizes, 
but seem to converge vs a value that is noticeably higher 
that 1/8 for q > in the limit L — > oo. 



The non-equilibrium phase q > 1 

For q > 1 the model has no stationary state and there- 
fore becomes sensitive to the initial conditions. For ex- 
ample, if the substrate is initially covered with several 
layers of particles of one type, particles of the other type 
are unable to attach so that the interface grows at con- 
stant velocity, just as in the single-species case. However, 
starting with a flat interface at zero height, the situation 
is totally different. Here one observes a slow coarsen- 
ing of competing three-dimensional islands of particles of 
both types. 

As for q > 1 the deposition rate is higher than the 
evaporation rate islands of both types tend to grow. The 
RSOS constraint, however, does not allow slopes greater 
than one and requires a line of unoccupied sites in be- 
tween islands of different types. The islands therefore 
have a pyramidal shape. There is a surface-tension-like 
effect, that makes large islands grow at the expense of 
small ones, but the coarsening in this phase is extremely 
slow. This is due to the fact that the borders of the is- 
lands are almost completely immobilized, as the particles 
in the bottom layer can evaporate only, if beforehand all 
the particles along the slope of the pyramidal island have 
evaporated. 

To quantify the slowdown of the dynamics we per- 
formed finite-size simulations with periodic boundary 
conditions. In these simulations, because of the finite 
system size, the coarsening process eventually displaces 
one of the particle types, allowing the interface to detach 
from the bottom layer and to propagate at constant ve- 
locity. As shown in Fig. Q21 the mean detachment time 
Td depends on the rates q$ and q and the system size L. 
The results suggest that the detachment time Td grows 
exponentially with the system size L, or even faster. 

Moreover, the simulations confirm the existence of two 
different regions (C) and (C) in Fig. [8] for q > 1, as 
already observed when measuring the Binder cumulant. 
In region (C) the formation of islands and the coarsen- 
ing process start immediately, while in region (C) one 
observes a low density of small short-lived islands for a 
long time until one or several large islands, that were 
generated by fluctuations, begin to grow. 

The different phenomenological properties in these re- 
gions can be explained as follows. For q > 1 the clusters 
grow by an influx of particles proportional to their area. 
If go < lj however, there is also a competing process: 
the tendency of particles to evaporate from the bottom 
layer. As such particles are preferentially located at the 
border of island the strength of this effect is proportional 
to the perimeter of an island. If qo is small enough the 
balance of these two processes defines a critical droplet 
size from where on islands start to grow. Although the 
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FIG. 12: (color online) Dependency of the displacement time 
Td on the rate qo and the system size L for q = 1.5. 



regions (C) and (C) seem to be well separated, we do 
not see evidence for a phase transition. 

For small system sizes L the qualitative different be- 
havior in these two regions even leads to quantitative in- 
fluence on the detach time (see Fig. [T2|). For high values 
of qo, i.e. in region (C), the asymptotically exponen- 
tial increase of the detachment time sets in already for 
small L. For small qo in region (C) the curves increase 
much faster for small L, leading to detachment times 
that are several orders of magnitude larger. Surprisingly 
the curves reach a local maximum followed by a decline 
and an eventual crossover to the asymptotic exponential 
increase. 

This behavior may be explained as follows: If the sys- 
tem size is smaller than the critical droplet size even sub- 
critical droplets can make the interface detach due to 
the periodic boundary conditions. The local maximum 
is reached when the critical droplet size approximately 
equals the system size. Increasing the system size fur- 
ther, the probability to generate a critical droplet grows, 
thereby reducing the detachment time. Here the dynam- 
ics is still governed by a single droplet which exceeds 
the critical size and then grows rapidly until the inter- 
face detaches. In much larger systems, however, it is 
likely that two or more supercritical droplets of different 
types are created before the first can spread out over the 
whole substrate. Then the limiting effect is the compe- 
tition between the pyramidal islands which explains the 
asymptotic exponential increase. 



the mapping of the model to a growth process with posi- 
tive and negative heights (see Fig. [3]) it is easy to see that 
in this case the potential p^|) is translationally invariant 
in the heights so that the dynamics of the model is equiv- 
alent to a freely evolving RSOS interface with average 
velocity zero without the hardcore wall at h = 0. At this 
point the surface does no more approach an equilibrium 
configuration, but rather it continues to roughen with a 
width increasing as ^/t. In fact, the density of unoccu- 
pied sites, that in this case is equivalent to the density 
of the roughening RSOS interface cut horizontally in the 
middle, decays like 



(27) 



Moreover, in a finite system starting at zero height the 
occupation balance is found to perform a random walk, 
hence 



0.500(5) 



for (d= 2,q = q = 1). (28) 



For qo > 1 translational invariance in the heights is bro- 
ken by a repulsive potential well at zero height. In the 
special case q = 1 and qo — > oo starting with a mono- 
layer of one type of particles, particles of the other type 
are unable to attach. The bottom layer acts as an effec- 
tive inert substrate on which particles of one type form 
clusters in the same way as they do in the RSOS model. 
As time proceeds, the interface roughens according to 
the Edwards- Wilkinson universality class. Therefore the 
width is expected to increase as 



t 



1/4 



(29) 



This is consistent with our numerical observation that 



sgn(6)) ~ t 



0.25(3) 



for [d = 2, q = 1, qo — * oo). 

(30) 



For go < 1 on the other hand one has an attractive 
potential barrier at zero height. When the line q = 1 
is approached from above the tendency of islands to 
grow diminishes as q decreases, and therefore the crit- 
ical droplet size gets larger and finally diverges as q = 1 
is approached. Starting with an empty lattice the model 
behaves similar as in the unordered phase (A). Even if 
the system is started with a few monolayers of particles 
fluctuations soon create holes in the deposit that grow 
until the substrate is almost completely exposed and the 
system again ends up in a situation that is similar to the 
unordered phase. 



CONCLUSIONS 

The line q = 1 

In the present work we have investigated how the phase 
Let us first consider the upper terminal point of the transition in a model for nonequilibrium wetting is af- 
phase transition line located at qo = q = 1. Using again fected by introducing two competing adsorbates with 



11 



identical properties. In this generalized model we have 
assumed that the two species of particles repel each other 
strongly. This leads to the formation of mutually re- 
pelling droplets, each consisting of only one type of par- 
ticles. 

In two dimensions we have identified an additional line 
of second-order phase transitions in the bound phase (see 
Fig. 4). It separates an ordered phase, where one of 
the two particle species takes over, from a disordered 
phase characterized by many small islands of different 
type. Defining appropriate order parameters and cumu- 
lants we have studied this transition in detail by numer- 
ical simulations. The origin of the phase transition can 
be explained in part by analyzing the partition sum in 
those parts of the phase diagram where detailed balance 
is valid. 

In the limit q — > 0, where the wetting process is gov- 
erned by a monolayer, this transition exhibits the same 
critical properties as the kinetic Ising model with heat 
bath dynamics. This is reasonable since in this case this 
islands are flat and their interior noise caused by evap- 
oration and deposition competes with ordering influence 
of surface tension at their boundaries. 

For < q < 1 higher layers are involved as well, 
giving the islands a three-dimensional shape. Never- 
theless, their thickness should remain finite and of the 
order ~l/hi(q) so that one would expect them to be- 
have asymptotically in the same way as in the monolayer 
case, suggesting that the entire curved line of transitions 
should belong to the Ising universality class. Contrarily, 
our numerical simulations seem to give evidence for con- 
tinuously varying exponents. The question whether these 
varying exponents are genuine or caused by crossover ef- 
fects is still open. A summary of the critical exponents 
can be found in Table HH 

As in the single-species case, the model undergoes a 
transition from a bound to an unbound phase at the 
threshold q = 1 (the horizontal line in Fig. [4|). How- 
ever, the critical properties for a system with an initially 
flat interface along this line are found to be different. In 
particular, there is the special point q = qo = 1, at which 
the two transition lines meet. It separates the line q = 1 
into two segments. For q$ > 1 the wetting transition is 
continuous, while for qo < 1 it is discontinuous. 

The unbound phase can be divided into two differ- 
ent regimes of different coarsening with an approximate 
borderline which seems to prolongate the curved transi- 
tion line from the bound phase. Nevertheless, we do not 
find evidence for a sharp phase transition between these 
regimes. In both regions the competition of two species 
slows down the dynamics. 

The study presented here has been primarily of theo- 
retical interest and was motivated by the question how 



an additional symmetry between competing adsorbates 
affects the critical properties of a wetting transition. As 
soon a this symmetry is broken one expects the process to 
cross over to an effective single-species behavior after suf- 
ficiently long time. It would be interesting to study the 
case p 1, where the process is out of equilibrium even 
in the bound phase. In this case the additional transition 
would still exist but probably it will no longer belong to 
the Ising universality class in the limit of small q. 
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